First, we should transfer this ODE into a difference equation form:
$$ x(t+\Delta t)-x(t)=-U_x'(x(t))\Delta t+\Delta t\xi $$where $\Delta t\xi$ denotes $\int_0^{\Delta t}\xi(s)ds$
Since $\xi(t)$ is a markov progress and the autocorrelation $\left<\xi(t_1)\xi(t_2)\right>=2\delta(t_1-t_2)$, so $\Delta t\xi \sim \mathcal N (0,2\Delta t)$, then we program it.
First, load denpendency and configure.
# mute logging in jupyter notebook
import Logging
Logging.disable_logging(Logging.Warn)
# load dependency
using Distributions, Plots, StatsPlots
Program the function randwalk for random walk simulation.
function randwalk(X₀::Vector{<:Real},T::Real;σ²=2,dUₓ=x->0,Δt=0.1)
Δtξ = Normal(0,√(σ²*Δt))
N = length(Xâ‚€)
ts = 0:Δt:T
M = length(ts)
X = zeros(N,M)
for (i,t) in enumerate(ts)
X[:,i] = t == 0 ? X₀ : X[:,i-1] .+ rand(Δtξ,N) .- dUₓ.(X[:,i-1]).*Δt
end
return (ts=ts,X=X)
end
function distanim(ts::AbstractRange,X::Matrix{<:Real};F=10,xbins=:fd,ylim=(0,0.5))
@gif for (f,(t,Xₜ)) in zip(ts,eachcol(X)) |> enumerate
(f-1) % F != 0 && continue
histogram(Xₜ,normalize=:pdf,bins=xbins,label=false)
title!("Particles distribution at t=$t")
ylims!(ylim)
xlabel!("\$X(t=$t)\$")
ylabel!("Probability Density")
end
end
For $U(x)=0$, the Brownian particle undergoes a random walk. Please carry out the simulation to obtain the particle trajectory and compute $\left<[x(t)-x(0)]^2\right>$, which is expected to be linear in t. Determine the diffusivity from the numerical data for $\left<[x(t)-x(0)]^2\right>$ and compare it with the theoretical value.
For $U(x)=0$, the equation of motion turns to this:
$$ \frac{dx}{dt}=\xi $$Solve this ode, we get:
$$ x(t)=x(0)+\int_0^t\xi(s)ds $$then
$$ \left<[x(t)-x(0)]^2\right>=2t $$thus the theoretical value of diffusivity $D=\frac12 \frac{d}{dt}\left<[x(t)-x(0)]^2\right>=1$.
Let number of particles $N=3000$, time duration $T=500$s, time step $\Delta t=0.1$s.
For simplicty, let $x(0)=0$ for all particle.
ts,X=randwalk(zeros(3000),500;σ²=2);
plotly()
plot(ts,X[end,:],label=false)
for i in 1:50
plot!(ts,X[i,:],label=false)
end
title!("Particle trajectory")
xlabel!("t")
ylabel!("X")